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A pseudo-Newtonian Hill problem based on a potential proposed by Artemova et al. [As- 
troph. J. 461 (1996) 565] is presented. This potential reproduces some of the general rela- 
tivistic effects due to the spin angular momentum of the bodies, like the dragging of inertial 
frames. Poincare maps, Lyapunov exponents and fractal escape techniques are employed to 
study the stability of bounded and unbounded orbits for different spins of the central body. 
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1 Introduction 



The Hill Problem was first formulated by Hill [1], in order to study the Moon- 
Earth-Sun system. This is a special case of the circular, planar restricted three- 
body problem, as described by Murray and Dermott [2] and Arnold [3], where the 
movement of the Moon around the Earth was just perturbed by a distant Sun. The 
Hill Problem is still applied in solar system models where bodies in nearly circular 
orbits are perturbed by other far away massive bodies, and is very useful in the 
study of the stellar dynamics. In many systems the Hill problem can be taken as a 
first approximation and can easily accommodate necessary modifications (see for 
instance Heggie [4]). The interaction of a Keplerian binary system with a normally 
incident circularly polarized gravitational wave can be represented by a Hill system, 
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as shown by Chicone et al. [5]. The Hill problem was proved to be non-integrable 
by Meletlidou et al. [6], and is chaotic, as shown by Simo and Stuchi [7]. 



The Hill Problem was first formulated in the context of the Newtonian dynamics. 
In this dynamics particles or spheres with or without spin have the same gravita- 
tional potential. However, on the General Relativity rotation modifies the metric 
and create effects like the dragging of inertial frames, or Lense-Thirring effect [8]. 
In General Relativity the exterior field of the simplest meaningful rotating body is 
described by the Kerr metric [9]. The Schwarzschild metric being the special case 
of the Kerr metric for spinless bodies. 

On a previous work [10], the Hill Problem is study in the framework of a pseudo- 
Newtonian potential that mimics some properties of the Schwarzschild metric, the 
Paczyhski-Wiita potential [11] 



GM 

O pw = , (1) 

r-r„ 



where G is the gravitational constant, M is the mass of the central body and r g = 
GM/c 2 is the Schwarzschild or gravitational radius (c is the speed of light). The 
aim of this work is to study the Hill problem taking into account the spin angular 
momentum of the central body, with a potential that mimics the dynamics of the 
Kerr metric. There are several potentials, originally used to study accretion disks 
around black holes, that fulfill this requirement. Some of this potentials are the 
Smerak-Karas potential [12] and the Mukhopadhyay potential [13]. However, the 
potential proposed by Artemova, Bjornsson and Novikov (ABN potential) [14] was 
chosen due to (i) its agreement with the last stable and marginally bound orbits, 
obtained from the Kerr metric itself; (ii) its simple form, a natural extension of the 
Paczyhski-Wiita potential. 

This work can be considered as a "zeroth order" approach to a general relativistic 
Hill problem with a spinning central body. Due to the complexity of the Einstein 
equations and its rigorous approximations, like the post-Newtonian expansions, it 
is worth to begin with this simplified approach to have an idea of the size of the 
quantities involved. These results can serve as a starting point for a more complete 
treatment of the problem. 

We shall compare the stability of orbits of the third body in the pseudo-Newtonian 
general relativistic simulation for different spins of the central body with parameters 
that are typical for a system formed by a supercluster, a galaxy and a star. In this 
system the influence of the spinning on the dynamics can be easily seen. 
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2 The ABN Potential 



Artemova, Bjornsson and Novikov [14] proposed two potentials to use in modeling 
accretion disks around a rotating black hole. They demanded that their potentials 
should have three properties: (i) the free-fall acceleration should have the form 
analogous to the Paczyhski-Wiita potential (equation (1)); (ii) the free-fall acceler- 
ation must tend to infinity when tends to the event horizon of the black hole; and 
(iii) the position of the extremum of the boundary condition function must coincide 
with the position of the last stable circular orbit in the exact relativistic problem. 
The boundary condition function reads 

fir) = 1 - T^r, (2) 
l{r) 

where l{r) = r 2 Q.(r) is the Keplerian specific angular momentum at radius r ( Q.(r) 
is the angular velocity being f2(r) = (GM/r 2 ) l/2 for the Newtonian gravitational 
potential). The constant /,„ is given by /,„ = l(r in ), where r in is the radius of the last 
stable orbit [14]. The first of the two the free-fall accelerations obtained is of the 
form 



GM 



ABN — 2 - fir V?' 

r l p{r — r\f 
or, in the potential form, 



(3) 
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(4) 



where r\ is the position of the event horizon. This position is determined from the 
angular momentum a by the exact expression from general relativity, 



ri = [l + (l-a 2 ) 1/2 ]r g . (5) 



The value of is given from the following equation: 

0= — -l (6) 



Where r in is again the last stable orbit. The exact position of this orbit is obtained 
from the general relativity (see Novikov and Frolov [15], equation 4.5.12), 

r in = [3 + Z 2 + [(3 - ZO(3 + Z l + 2Z 2 )] 1/2 ]r„ (7) 
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Table 1 

Values of r m b / r g for the ABN potential and Kerr geometry for various values of a 



a 





0.1 


0.3 


0.5 


0.7 


0.998 


$>ABN 


4 


3.797 


3.370 


2.904 


2.376 


1.069 


Kerr 


4 


3.797 


3.373 


2.914 


2.395 


1.091 


a 





-0.1 


-0.3 


-0.5 


-0.7 


-0.998 


®ABN 


4 


4.198 


4.578 


4.942 


5.288 


5.746 


Kerr 


4 


4.198 


4.580 


4.949 


5.308 


5.825 



Z, = 1 + (1 - a 2 ) 1/3 [(l + a) 1 ' 3 + (1 - a) 1/3 ], (8) 

Z 2 = (3a 2 + Z 2 ) 1/2 (9) 

where the upper and the lower signs of the equation (7) are for co-rotation and 
counter-rotation respectivel>Q- Note that the parameters r in and ft depend only on 
the angular momentum a. 

The radius of the last stable orbit, r s matches exactly with the last stable orbit in 
Kerr geometry, as demanded on the construction of this potential. In Table 1 the 
values of the marginally bound orbit, r mb are listed for the potential above and the 
Kerr geometry for various values of a. Note that they are in a good agreement 
with each other. Mukhopadhyay [13] has pointed that using the ABN potential for 
negative values the error in r b may be upto 500%. However, if the correct equation 
for counter-rotation is used (equation (7) with the lower signal), the error in r b is 
less than 2%. 

For the second potential a different expression for the free-fall acceleration F ABN is 
obtained. We shall use (4) because it is singular in the event horizon r = r\, like the 
Paczyhski-Wiita potential. The second potential proposed by ABN does not have 
the above mentioned property. 



3 Modified Hill Problem 



The Hill problem is a special case of the circular, planar restricted three-body prob- 
lem, as mentioned before. In this problem there is a system of two massive bodies 
with masses nt\ and m 2 , in circular orbits around their center of mass and a third 
massless body moving under influence of this system without perturbing it. In the 
Hill problem the body with mass m\ is such that mi » m 2 and is far away of the 
system, so it constitutes just a perturbation for the two-body system formed by the 

1 In Ref. [14] equation (7) appears only with the minus sign, but the correct form is the 
above, as in Novikov and Frolov [15]. 
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Moon 




1-U 

Fig. 1. The planar, circular restricted three-body problem, on a inertial coordinate system 
(£, tj) and on a rotating coordinate system (x, y). Note that on the rotating coordinate system 
the positions of the two massive bodies are fixed. 

body with mass m 2 and the massless body. We can choose units of mass such that 
G(nii +m 2 ) = 1. In this way we can take the units of mass of the two massive bodies 
respectively as 1 -fj. and /i, where [x = Gm 2 . We take units of distance such that the 
distance between the two massive bodies is R = 1 and units of time such that the 
angular velocity of the rotating frame in which the two massive bodies are fixed is 
a> = 1 (see Fig. 1). 

After placing the origin of the coordinate system in the body of mass // and con- 
sidering the motion only in a disk of radius // 1/3 , we obtain the Newtonian Hill 
equations [7], 



x = 2y + 3x — -, 
y = -2x- — , 



(10) 

(11) 



where r = sjx 2 +y 2 . Replacing the Newtonian potential by the ABN potential we 
obtain the modified Hill equations, 



x = 2y + 3x - 



x 



(12) 
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These equations can also be written as 



x - 2 y = - 
y + 2x = - 



dU 



ABN 



dx 
9Uabn 



By ' 

where U AB n is the modified Hill potential 



(14) 
(15) 



3 2 1 

Uabn - ~~z x ~ 77, — 7~T* \ 

2 (J3- l)r* [\r- r\ 



- 1 



3 2 ^ 
-^X 2 + d>ABN, 



(16) 



and r* = [i 1/3 ri, in units such that the separation between the two massive bodies 
is R = 1 , like in the Newtonian case. 



The modified Jacobi constant, in this case, is given by 



C 



J ABN 



1 / -2 -2\ 3 2 1 



r - r: 



- 1 



(17) 



4 Stability of orbits 



We shall use for the parameter r* the value r* = 5.10 6 , that is a typical value for 
a system formed by a supercluster, a galaxy and a star. We compare the modified 
Newtonian systems with the angular momentum a varying from -0.5 to 0.5. The 
values of f3 vary, respectively, from 3.05 to 1.27. 



4.1 Fixed points analysis 



The Newtonian Hill problem has two well-known fixed points, the Lagrangian 
points L\ and L 2 . These points are of saddle-center type, and are located at the po- 
sitions (+-^1/3, 0). For the ABN potential it can be guessed that new fixed points, 
in particular saddle points, arise due to the r - r\ dependence of the denominator. 
However, if the parameter r\ is small, the nature of the fixed points remains un- 
changed. The equation for the x component of the fixed points (obtained by setting 
x, y, 'x, y equal to zero in equations (12) and (13) and replacing y = 0) is, 



\xt P (\x\-r\f= 1 -. 



(18) 



If r\ is small compared to x (in this work r\ is of order 10 6 while x scales as unity) 
the equation (18) reads, approximately, 



'4-4)4 



(19) 



For the values of Cj so that the system is bounded, the value of the x component 
of the fixed points is also bounded. There are one real and two complex conjugate 
solutions for this approximate equation. For higher values of r\ the exact equation 
18) can have more than one solution, they may possibly lead to interesting results, 
but these systems may lack of physical significance. 

The Jacobian of the system calculated at the fixed points reads 



J F (x f ,y f ) 



10 
1 

9+3 ^ 2 

-3-2 



(20) 



The associated characteristic polynomial is 



(* + 3)(*-9-*5fy 



o. 



(21) 



The eigenvalues are 



A = +/V3 



A=±J9+3fi 



\*f\ ~ r \ 



(22) 
(23) 



For r* small enough [ r\ < 3\xf\/(3 +/?)] there are two eigenvalues purely imaginary 
and two real, so the fixed points are of type saddle-center, just as in the Newtonian 
Hill problem. 



4.2 Poincare sections 



The orbits in the Newtonian as well in the modified Newtonian Hill problem are the 
solutions of a four-dimensional dynamical system with variables (x,y,x, y). Since 
we have an integral of motion, Cjabn, the motion is reduced to a three-dimensional 
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system, we can take (x, x,y) as independent variables.We shall study surface of 
section (Poincare sections) evaluating the orbits for different values of the Jacobi 
constant and registering the crossings of the hypersurface y = with y > 0. 

The results for Cj = -2.17 (typical value for a bounded system) are shown on the 
Fig. 2. In Fig. 3 there is a magnification of the right portion for the values a = 0.5 
and a = -0.5. We see that, as the angular momentum a increases in modulus (from 
(a) to (c) and from (d) to (f) in Fig. 2) some Kolmogoro v- Arnold-Mo ser (KAM) 
tori are deformed and others destroyed, indicating the transition of the system from 
regular to chaotic behaviour. This transition is more evident for the sequence of 
negative values of a, i.e., the sequence from (d) to (f), so the destruction of KAM 
tori is faster for counter-rotation. It means that, for the counter-rotation case, a 
larger region of the phase space is chaotic when compared with the associated co- 
rotation system (same angular velocity a in modulus). It does not means necessarily 
that the counter-rotating systems are more unstable than the co-rotating ones, since 
the orbits are still bounded by KAM tori that are not destructed by the perturbation. 
Besides, the dependence on initial conditions, the main characteristic of chaotic 
systems, still must be analyzed by appropriate tools, like Lyapunov exponents. This 
analysis can decide if a system is more chaotic than other [16]. 



4.3 Lyapunov exponents 



To analyze quantitatively the orbits stability we study the Lyapunov exponents for 
the systems above described. The Lyapunov characteristic number (A) is defined as 
the double limit, 



A = lim 

5o->0 



log05/<5 ) 



(24) 



where So and 8 are the deviation of two nearby orbits at times and t respectively 
(see Alligood et al. [17]). We get the largest A using the technique suggested by 
Benettin et al. [18] and the algorithm of Wolf et al. [19]. 

The Lyapunov exponents are not absolute, but dependent on the choice of the time 
scale. We recall that we have fixed the time scale by the requirement that a> = 1 . 
This defines a time unit that is natural to each particular system at it is given in 
terms of the characteristic period. In this work the analysis is made by varying only 
the angular momentum a, so the direct comparison between the different Lyapunov 
exponents is valid. Each coefficient was computed until convergence is reached. 
To achieve this precision the system of equations was integrated for at least one 
hundred thousand periods, as shown in Fig. 4. The initial conditions used to perform 
the integration must be chosen in the bounded region of the system. To estimate the 
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0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 



Fig. 2. Poincare sections for different values of the angular momentum a. Positive values 
a) a - 0.1 , b) a - 0.3 and c) a = 0.5 and negative values d) a = -0.1, e) a - -0.3 and 
f) a = -0.5. The rate of destruction of KAM tori is greater for the negative values. Note 
that, for larger values of a, new islands appear in the regular region, indicating transition to 
chaos. 

limit of this region we use the Lagrangian points of the Newtonian Hill system, 
given numerically by (+0.69,0). For safety we chose the initial position (0.3,0). 
The initial velocity is obtained from the Jacobi constant Cj = -2.17, choosing 
x = 0. As can be seen in the subsection 4.1 the eigenvalues obtained from the 
Jacobian of the system remain small, so the system is non-stiff. The calculations 
are performed with the aid of the Burlisch-Stoer method with step control, that 
works well for non-stiff systems, and the error due to the integration is proportional 
to the tolerance imposed (10 10 ). This relative error is kept small enough, so the 
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Fig. 3. Magnifications of areas in the Poincare sections, at left for angular momentum 
a - 0.5 and at right for angular momentum a = -0.5. Note that for a = -0.5 the tori are 
more affected by the perturbation. 

errors of the coefficients are basically estimate from the fluctuations at the end of 
each evolution (Fig. 4). The absolute error, according to these fluctuations, is on 
order of 10 8 . 

The Lyapunov exponents obtained for this system are shown in Fig. 5. As can be 
seen from this figure, negatives value of a are associated to a larger Lyapunov ex- 
ponent when compared with their correspondent positive values. It means that two 
neighbour orbits separate from each other faster in the counter-rotating system than 
in the co-rotating system, meaning that the counter-rotating system is more unsta- 
ble. Nevertheless, this dependence on parameter a is very small for real systems. 
The reason is the fact that the apparent event horizon is very small when compared 
with typical distances in the system, and the angular momentum a has only little 
influence on this horizon as already pointed by Bardeen, Press and Teukolsky [20]. 
We shall study this dependence in a future work taking another model for the Hill 
problem closer to the exact general relativistic system. 

d) Fractal escape and fractal dimension 

The Poincare sections and Lyapunov exponents were obtained for values of the Ja- 
cobi constant such that the systems are bounded. For values larger than C Jbounded = 
-2.16 the systems are unbounded and the third body can escape by two different 
routes [10]. For open systems that have more than one route to escape we can apply 
the fractal escape technique used by Moura and Letelier [21] in the study of the 
classical Henon-Heiles problem. In this method the basins of the escape routes are 
obtained for a set of initial conditions. For chaotic systems, we have the existence 
of fractal basin boundaries (FBB) indicating a great instability of the orbits. In our 
case we chose a subset of the accessible phase space at a fixed Jacobi constant, de- 
fined by a segment \x\ < a,y = b and < 6 < 2n, where a and b are constants to be 
chosen appropriately and 8 is the angle that defines the direction of the velocity with 
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100 1000 10000 le+05 le+06 

Integration time (t) 

Fig. 4. Time evolution of the separation exponent. It converges to the corresponding Lya- 
punov characteristic number (A) for large t. 




0.14082 



Angular Momentum (a) 



Fig. 5. Lyapunov exponents for different values of a. 
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respect the jc-axis. Then the trajectories are integrated numerically, we have three 
different cases: (i) the body escapes to x —> +00, (ii) it escapes to x —> -00, and (iii) 
the particle does not scape during the integration time. We take the integration time 
long enough to assure that our results be consistent. 

To show the difference between systems with different values of the parameter a, we 
calculate the dimension of the basin boundaries obtained for the different systems. 
The dimension used is the box-counting dimension that can be easily obtained, see 
for instance, Ott [22] and Grebogi et al. [23]. If we displace a determined point of a 
basin to another on a distance e, the probability that this new initial condition does 
not belong to the same basin of the old one is, for small e, P(e) oc e D ~ d , where D 
is the dimension of the set (2 in our case) and d is the box counting dimension, 
also called exterior or fractal dimension when not integer. In order to calculate this 
fractal dimension, for several values of e, we displace the x coordinate of all the 
points from one of the basins (x —> +00), and count the number of points that does 
not belong to the same basin. Then we compute the the fraction of numbers that 
does not belong to the same basin, P(e). We plot ln.P(e) in function of In 6. The 
inclination of the straight line gives us D - d. 

The values for the fractal dimension obtained are shown in the Fig. 6, with error 
of order 3.10 3 . This uncertainty is mainly due to the error in the computation of 
the line slope. It can be seen that the fractal dimension for negative values of a are 
larger than the correspondent positive values. The fractal dimension for unbounded 
systems have a close relationship with the chaoticity of the systems [10], [21], larger 
fractal dimensions are related to systems that are more unstable. In this sense we 
can conclude that counter-rotating systems are more unstable than their correspon- 
dent co-rotating cases. The validity of this result is limited, as it is applicable only 
to unbounded systems and the error in calculations are of the order of the variation 
of data. For bounded systems the analysis of the Lyapunov exponents in the previ- 
ous section confirms the validity of this result. The values of d are computed only 
for a few values of a. Is possible that this small density of points can hide a more 
complex behaviour for this system. Unfortunately to achieve the same density of 
points used for the Lyapunov exponents is still prohibitive due to the time taken to 
perform the calculations. We shall improve this method in our future works. 



5 Conclusions 

In this work we show that we can simulate the dragging of inertial frames by a 
pseudo-Newtonian potential. We also show, by Poincare maps, Lyapunov expo- 
nents and fractal escape techniques, that there is a dependence of the stability of 
the orbits on the spin angular momentum of the central body. The bounded and 
unbounded systems where the movement of particle around the central body is 
opposite to its spin (counter-rotating) are more unstable than systems where the 
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Fig. 6. Fractal dimensions for different values of a. 

two rotations are in the same direction (co-rotating). This preliminary result is in 
accord with previous studies of the stability of orbits of particles moving around 
spinning centers of attraction [24], [25]. This effect is small when compared with 
the influence of the position of the event horizon [10]. Otherwise, it can have a 
larger influence on the Lyapunov exponents. In a future work this feature will be 
studied in detail to compare the chaoticity of co-rotating and counter-rotating or- 
bits in a full relativistic system, as geodesies in a Kerr black hole with halos [24] or 
multipolar deformations [25]. 
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